WindGonzalezLongatt Subroutine

SUBROUTINE WindGonzalezLongatt &
!
(speed, dir, dem, grid, winddir)

USE Units, ONLY: &
!Imported parameters
degToRad, radToDeg, pi

IMPLICIT NONE

!Arguments with intent (in):
TYPE (ObservationalNetwork), INTENT(IN) :: speed !windspeed observations at stations 
TYPE (ObservationalNetwork), INTENT(IN) :: dir !wind direction observations at stations 
TYPE (grid_real), INTENT(IN) :: dem !digital elevation model

!Arguments with intent (inout):
TYPE (grid_real), INTENT (INOUT) :: grid !wind speed grid [m/s]
TYPE (grid_real), OPTIONAL, INTENT (INOUT) :: winddir !wind direction [deg]


!Local declarations
TYPE (ObservationalNetwork) :: temp_speed, temp_dir, active_speed, active_dir
TYPE (ObservationalNetwork) :: uwind, vwind
INTEGER :: count_speed, count_dir, i, j
TYPE (grid_real) :: uwind_grid !zonal wind [m/s]
TYPE (grid_real) :: vwind_grid !meridional wind [m/s] !OK
REAL (KIND = float) :: vxcalc !zonal wind with orograpcic correction applied [m/s]
REAL (KIND = float) :: vycalc !meridional wind with orograpcic correction applied [m/s]
REAL (KIND = float) :: corr1, corr2, cellsize
LOGICAL :: validNorth, validSouth, validEast, validWest


!----------------end of declarations-------------------------------------------

!check for available measurements. Both speed and direction must exist
!first check nodata. If only speed or direction are missing, set both to missing 
! assume dir and speed stations are in memory in the same order
CALL CopyNetwork (speed,temp_speed)
CALL CopyNetwork (dir,temp_dir)

DO i = 1, temp_speed %countObs
  IF (temp_speed % obs (i) % value == temp_speed % nodata) THEN
      !set dir to nodata
      temp_dir % obs (i) % value = temp_dir % nodata
  ELSE IF (temp_dir % obs (i) % value == temp_dir % nodata) THEN
      !set speed to nodata
      temp_speed % obs (i) % value = temp_speed % nodata
  END IF
END DO

!now remove nodata
CALL ActualObservations (temp_speed, count_speed, active_speed)
CALL ActualObservations (temp_dir, count_dir, active_dir)

!initialize zonal and meridional component network
CALL CopyNetwork (active_speed, uwind)
CALL CopyNetwork (active_speed, vwind)

!initialize grid
CALL NewGrid (uwind_grid, grid)
CALL NewGrid (vwind_grid, grid)

!compute u and v at stations
DO i = 1, active_speed % countObs
  uwind % obs (i) % value = - active_speed % obs (i) % value * &
                          SIN ( active_dir % obs (i) % value * degToRad)
  vwind % obs (i) % value = - active_speed % obs (i) % value * &
                          COS ( active_dir % obs (i) % value * degToRad)
END DO

!interpolate u and v grid 
CALL IDW (network=uwind, grid=uwind_grid, method=1, power=2.0, near=5)
!CALL BarnesOI (network=uwind, grid=uwind_grid, gammapar = 0.6)
     
CALL IDW (network=vwind, grid=vwind_grid, method=1, power=2.0, near=5)
!CALL BarnesOI (network=vwind, grid=vwind_grid, gammapar = 0.6)


!apply orographic correction
DO i = 1, grid % idim
  DO j = 1, grid % jdim
    IF (grid % mat (i,j) /= grid % nodata) THEN

       !set cellsize
       cellsize = (CellArea (grid,i,j)  ) ** 0.5 ![m]
       
       !check for neighbour cells. Nodata and out
       !of grid cells must not be considered
       !for orographic correction     
       validNorth = CellIsValid (i-1, j, grid)
       validSouth = CellIsValid (i+1, j, grid)
       validEast  = CellIsValid (i, j+1, grid)
       validWest  = CellIsValid (i, j-1, grid)

       !correct zonal (along x) component
       IF (validEast .AND. validWest)  THEN !apply correction where east and west data are defined

          !west correction factor component
          corr1 = (dem % mat (i,j)  -  dem % mat (i,j-1)) * &
                (uwind_grid % mat (i,j-1) - uwind_grid % mat (i,j)  ) / cellsize
          
          !east correction factor component
          corr2 = (dem % mat (i,j+1)  -  dem % mat (i,j)) * &
                (uwind_grid % mat (i,j) - uwind_grid % mat (i,j+1)  ) / cellsize
          vxcalc = uwind_grid % mat(i,j) + 0.5 * ( corr1 + corr2)
       ELSE !do not apply correction
          vxcalc = uwind_grid % mat(i,j)
       END IF

       !correct meridional (along y) component
       IF (validNorth .AND. validSouth)  THEN  !apply correction north and south data are defined

          !south correction factor component
          corr1 = (dem % mat (i,j)  -  dem % mat (i+1,j)) * &
                (vwind_grid % mat (i+1,j) - vwind_grid % mat (i,j)  ) / cellsize
          
          !north correction factor component
          corr2 = (dem % mat (i-1,j)  -  dem % mat (i,j)) * &
                (vwind_grid % mat (i,j) - vwind_grid % mat (i-1,j)  ) / cellsize
          vycalc = vwind_grid % mat(i,j) + 0.5 * ( corr1 + corr2)

       ELSE !do not apply correction
          vycalc = vwind_grid % mat(i,j)
       END IF  
       
       !compute wind direction if required
       IF (PRESENT (winddir)) THEN
           ! Some compilers do not allow both u and v to be 0.0 in
           ! the atan2 computation.
           IF ( ABS(vxcalc ) < 1e-10) vxcalc = 1e-10
           
           winddir % mat (i,j) = ATAN2(vxcalc,vycalc)
  
           IF (winddir % mat (i,j) >= pi) THEN
              winddir % mat (i,j) = winddir % mat (i,j) - pi
           ELSE
              winddir % mat (i,j) = winddir % mat (i,j) + pi
           END IF
           !convert radians to deg
           winddir % mat (i,j) = winddir % mat (i,j) * radToDeg
       ENDIF
           
        !windspeed
        grid % mat (i,j) = SQRT(vxcalc**2. + vycalc**2.)
        
    END IF
  END DO
END DO


!deallocate memory
CALL DestroyNetwork (temp_speed)
CALL DestroyNetwork (temp_dir)
CALL DestroyNetwork (active_speed)
CALL DestroyNetwork (active_dir)
CALL DestroyNetwork (uwind)
CALL DestroyNetwork (vwind)

CALL GridDestroy (uwind_grid)
CALL GridDestroy (vwind_grid)

RETURN
END SUBROUTINE WindGonzalezLongatt